######################################################################
# load data
rm(list = ls())
setwd("~/Dropbox/VMR_Project")
load("./data/raw_all.RData")

# exclude repeated rows in the raw data and export to a new list raw2

raw2 = list()

raw2$d20160818_2 = `20160818_2.XLS`[1:(1094590/2)*2, -5]
raw2$d20160818_2 = raw2$d20160818_2[1:547200,]

raw2$d20160818_4 = `20160818_4.XLS`[1:(1094590/2)*2, -5]
raw2$d20160818_4 = raw2$d20160818_4[1:547200,]

raw2$d20160819_4 = `20160819_4.XLS`[1:(1094590/2)*2, -5]
raw2$d20160819_4 = raw2$d20160819_4[1:547200,]

raw2$d20160819_6 = `20160819_6.XLS`[1:(1094400/2)*2, -5]
raw2$d20160819_6 = raw2$d20160819_4[1:547200,]
    
raw2$d20160922_1 = `20160922_1.XLS`[1:(1094590/2)*2, -5]
raw2$d20160922_1 = raw2$d20160922_1[1:547200,]

raw2$d20160922_2 = `20160922_2.XLS`[1:(1094590/2)*2, -5]
raw2$d20160922_2 = raw2$d20160922_2[1:547200,]

raw2$d20160926_2 = `20160926_2.XLS`[1:(1094592/2)*2, -5]
raw2$d20160926_2 = raw2$d20160926_2[1:547200,]

raw2$d20160926_3 = `20160926_3.XLS`[1:(1094590/2)*2, -5]
raw2$d20160926_3 = raw2$d20160926_3[1:547200,]

raw2$d20160927_1 = `20160927_1.XLS`[1:(1094598/2)*2, -5]
raw2$d20160927_1 = raw2$d20160927_1[1:547200,]

raw2$d20160927_2 = `20160927_2.XLS`[1:(1094598/2)*2, -5]
raw2$d20160927_2 = raw2$d20160927_2[1:547200,]

raw2$d20160927_3 = `20160927_3.XLS`[1:(1094590/2)*2, -5]
raw2$d20160927_3 = raw2$d20160927_3[1:547200,]

raw2$d20160927_4 = `20160927_4.XLS`[1:(1094590/2)*2, -5]
raw2$d20160927_4 = raw2$d20160927_4[1:547200,]

raw2$d20161004_1 = `20161004_1.XLS`[1:(1094590/2)*2, -5]
raw2$d20161004_1 = raw2$d20161004_1[1:547200,]

raw2$d20161004_2 = `20161004_2.XLS`[1:(1094590/2)*2, -5]
raw2$d20161004_2 = raw2$d20161004_2[1:547200,]

raw2$d20161004_3 = `20161004_3.XLS`[1:(1094590/2)*2, -5]
raw2$d20161004_3 = raw2$d20161004_3[1:547200,]

raw2$d20161004_4 = `20161004_4.XLS`[1:(1094590/2)*2, -5]
raw2$d20161004_4 = raw2$d20161004_4[1:547200,]

raw2$d20161004_5 = `20161004_5.XLS`[1:(1094590/2)*2, -5]
raw2$d20161004_5 = raw2$d20161004_5[1:547200,]

raw2$d20161004_6 = `20161004_6.XLS`[1:(1094590/2)*2, -5]
raw2$d20161004_6 = raw2$d20161004_6[1:547200,]


# export 30 secs before and 30 secs after light-on to a new list, current3
current3 = list()
for (ij in 1:18) {
    current3[[ij]] = raw2[[ij]][(1+96*1769):(1830*96),c(4,6,13,16,19)]
    current3[[ij]]$mean = apply(current3[[ij]][,3:5], 1, sum)
    cat(ij, " ")
}
names(current3) = names(raw2)
save(current3, file = './data/lightOn30.RData')


# export 60 secs before and 120 secs after light-off to a new list, current2

current2 = list()
for (ij in 1:18) {
    current2[[ij]] = raw2[[ij]][(1+96*5340):(5520*96),c(4,6,13,16,19)]
    current2[[ij]]$mean = apply(current2[[ij]][,3:5], 1, sum)
    cat(ij, " ")
}
names(current2) = names(raw2)

rm(raw2,ij)
rm(list = ls()[1:18])
rm(current3)

save.image("./data/current_all.RData")


######################################################################
# load data
rm(list = ls())
load("./data/current_all.RData")
current2.avg = current2$d20160818_2[,3:6]
for(i in 2:18){
    current2.avg = current2.avg + current2[[i]][3:6]
}
current2.avg = current2.avg / 18
current2.avg$genotype = current2$d20160818_2$genotype
current2.avg$start = current2$d20160818_2$start

tmp2 = as.matrix(current2.avg[,1:4]) # convert to matrix
zbc = array(NA, dim = c(96, 4, 180), 
            dimnames = list(current2.avg$genotype[1:96], names(current2.avg)[1:4], 
                            5340:5519))
for (ij in 1:180) {
    zbc[,,ij] = tmp2[96*(ij-1)+1:96,]
    # cat(ij, ' ')
}
rm(tmp2, i, ij)


subtotals = sapply(1:180, function(k) t(aggregate(current2.avg[96*(k-1)+1:96,1:4], 
                                                  by=list(current2.avg$genotype[1:96]), 
                                                  mean)[,2:5]))
row.names(subtotals) = c(paste0('Q_', names(current2.avg)[1:4]), 
                         paste0('R_', names(current2.avg)[1:4]))
colnames(subtotals) = 5340:5519 - 5400

func1 = function(x) {sqrt( var(x) / 48 )}
subtotals.sd = sapply(1:180, function(k) t(aggregate(current2.avg[96*(k-1)+1:96,1:4], 
                                                  by=list(current2.avg$genotype[1:96]), 
                                                  func1)[,2:5]))
row.names(subtotals.sd) = c(paste0('Q_', names(current2.avg)[1:4]), 
                         paste0('R_', names(current2.avg)[1:4]))
colnames(subtotals.sd) = 5340:5519 - 5400

pic.show = list()
library(ggplot2)
for(ij in 1:4){
    df.tmp = data.frame(value = c(subtotals[ij,], subtotals[ij+4,]), 
                        se = c(subtotals.sd[ij,], subtotals.sd[ij+4,]),
                        sec = rep(5340:5519, 2),
                        genotype = rep(c("Q344X", "Rho"), each = 180))
    rib.off = aes(ymax = value + se, ymin = value - se)
    pic.show[[ij]] = ggplot(data=df.tmp, aes(x=sec, y=value, colour=genotype)) + 
        geom_line() + ylab("Distance Travelled (cm)") + xlab("Second (start)") +
        geom_ribbon(rib.off, fill = "grey70", alpha = I(15/100), linetype=0) +
        geom_vline(xintercept=5400, linetype="dashed") + theme_light() + 
        ggtitle(names(current2.avg)[ij]) + theme(plot.title = element_text(hjust = 0.5, size=16))
    
}

library(gridExtra)
grid.arrange(pic.show[[1]], pic.show[[2]], pic.show[[3]], pic.show[[4]], ncol = 2)

######################################################################
# plots of single excel 
rm(list = ls())
load("./data/current_all.RData")

current2.avg = current2$d20161004_4
subtotals = sapply(1:180, function(k) t(aggregate(current2.avg[96*(k-1)+1:96,3:6], 
                                                  by=list(current2.avg$genotype[1:96]), 
                                                  mean)[,2:5]))
row.names(subtotals) = c(paste0('Q_', names(current2.avg)[3:6]), 
                         paste0('R_', names(current2.avg)[3:6]))
colnames(subtotals) = 5340:5519

func1 = function(x) {sqrt( var(x) / 48 )}
subtotals.sd = sapply(1:180, function(k) t(aggregate(current2.avg[96*(k-1)+1:96,3:6], 
                                                     by=list(current2.avg$genotype[1:96]), 
                                                     func1)[,2:5]))
row.names(subtotals.sd) = c(paste0('Q_', names(current2.avg)[3:6]), 
                            paste0('R_', names(current2.avg)[3:6]))
colnames(subtotals.sd) = 5340:5519

pic.show = list()
library(ggplot2)
for(ij in 1:4){
    df.tmp = data.frame(value = c(subtotals[ij,], subtotals[ij+4,]), 
                        se = c(subtotals.sd[ij,], subtotals.sd[ij+4,]),
                        sec = rep(5340:5519, 2),
                        genotype = rep(c("Q344X", "Rho"), each = 180))
    rib.off = aes(ymax = value + se, ymin = value - se)
    pic.show[[ij]] = ggplot(data=df.tmp, aes(x=sec, y=value, colour=genotype)) + 
        geom_line() + ylab("Distance Travelled (cm)") + xlab("Second (start)") +
        geom_ribbon(rib.off, fill = "grey70", alpha = I(15/100), linetype=0) +
        geom_vline(xintercept=5400, linetype="dashed") + theme_light() + 
        ggtitle(names(current2.avg)[ij+2]) + 
        theme(plot.title = element_text(hjust = 0.5, size=16))
    
}
library(gridExtra)
grid.arrange(pic.show[[1]], pic.show[[2]], pic.show[[3]], pic.show[[4]], ncol = 2)

######################################################################
# load light data
setwd("~/Dropbox/VMR_Project")
par(mfrow=c(2,2))
light1 = read.table('./data/light/Well_A1.txt', header = F, sep = '\t', skip = 10)
plot(light1, type = 'l', main = "Well A1", 
     xlab = "Wavelength(nm)", ylab = "Irradiance(uW/cm^2/nm)")

light2 = read.table('./data/light/Well_E6.txt', header = F, sep = '\t', skip = 10)
plot(light2, type = 'l', main = "Well E6", 
     xlab = "Wavelength(nm)", ylab = "Irradiance(uW/cm^2/nm)")

light3 = read.table('./data/light/Well_D9.txt', header = F, sep = '\t', skip = 10)
plot(light3, type = 'l', main = "Well D9", 
     xlab = "Wavelength(nm)", ylab = "Irradiance(uW/cm^2/nm)")

light4 = read.table('./data/light/Well_G4.txt', header = F, sep = '\t', skip = 10)
plot(light4, type = 'l', main = "Well G4", 
     xlab = "Wavelength(nm)", ylab = "Irradiance(uW/cm^2/nm)")

######################################################################
# load light-on data
rm(list = ls())
setwd("~/Dropbox/VMR_Project")
load('./data/lightOn30.RData')

# current3.avg = current3[[1]][,3:6]
# 
# for(i in 2:18){
#     current3.avg = current3.avg + current3[[i]][3:6]
# }
# current3.avg = current3.avg / 18
# current3.avg$genotype = current3$d20160818_2$genotype
# current3.avg$start = current3$d20160818_2$start

zbc4 = array(NA, dim = c(96, 4, 61, 18), 
            dimnames = list(current3$d20160818_2$genotype[1:96], 
                            c('inadist', 'smldist', 'lardist', 'mean'), 
                            1769:1829, 1:18))


for(ii in 1:18){
    cat(ii, ' ')
    for(ij in 1:61){
        for(jj in 1:4){
            zbc4[,jj,ij,ii] = as.matrix( current3[[ii]][96*(ij-1)+1:96,2+jj] )
        }
    }
}


zbc4.30BeforeLightOn = zbc4[,,1:30,]

dim(zbc4.30BeforeLightOn)

# baseline: 30 secs before light onset, 4 variables x 18 batches
baseline = matrix(NA, nrow = 4, ncol = 18)

for(ij in 1:4){
    for(jj in 1:18){
        baseline[ii, jj] = apply()
    }
}
baseline = sapply(1:96, function(x) apply(zbc3.30BeforeLightOn[x,,,1],1,mean))
# baseline1 for all files except 20160927_1
baseline1 = baseline
# baseline2 for file 20160927_1
baseline2 = baseline1[, c(15:96,1:14)]
save(baseline1, baseline2, file = './data/baseline.RData')
